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In this work we derive an analytic expression for the Kolmogorov-Sinai entropy of dilute wet 
granular matter, valid for any spatial dimension. The grains are modelled as hard spheres and the 
influence of the wetting liquid is described according to the Capillary Model, in which dissipation is 
due to the hysteretic cohesion force of capillary bridges. The Kolmogorov-Sinai entropy is expanded 
in a series with respect to density. We find a rapid increase of the leading term when liquid is added. 
This demonstrates the sensitivity of the granular dynamics to humidity, and shows that the liquid 
significantly increases the chaoticity of the granular gas. 
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I. INTRODUCTION 

The field of granular physics has undergone consider- 
able progress in recent times fT, '2*1 . As part of soft mat- 
ter physics, granulates have inspired the development of 
non-equilibrium statistical mechanics Its poten- 

tial to the foundation of physics can hardly be over esti- 
mated, since granular gases provide a road away from the 
well-developed Boltzmann-Enskog theory of conservative 
gases towards dissipative systems far from thermal equi- 
librium. In connection with geophysics, some aspects of 
landslides may be understood in terms of a solid-liquid 
phase transitions of wet granular matter [^, and 
wet granular gases are of technological relevance in gran- 
ulators, pelletizers, and other instances in process engi- 
neering. 

Wet granular gases are systems consisting of meso- 
scopic particles and a liquid phase wetting the particles. 
Despite their importance, the theory of wet granular mat- 
ter is still nascent. There is a growing number of exper- 
imental [1] and numerical work Q on this subject, but 
the hysteretic nature of the liquid bridge interaction was 
not taken into account in the modelling. We stress that 
the attraction force mediated by capillary bridges is not a 
function of distance but depends on the collision history. 
The theory of wet granular matter advanced with recent 
simulation and models describing the free cooling state 
fldi [T]| . To the best of our knowledge, the hysteretic 
dissipative dynamics of wet granular matter was treated 
analytically first in [l^]- In this article we elaborate on 
this approach which treats the wet granulate as a com- 
plex dynamical system and uses powerful tools available 
in this area. Such is the Lyapunov spectrum. 



It gives the rate of exponential divergence or conver- 
gence of two equal copies of the system in phase space. 
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STj{t) = r^-^^(i) - Tf\t), with perturbed initial condi- 
tions (5rj(0). A positive Lyapunov exponent indicates 
chaotic behavior, i.e. sensitive dependence on the initial 
conditions [l^. Since we are dealing with a closed sys- 
tem the sum of all positive Lyapunov exponents equals 
the Kolmogorov-Sinai entropy (KSE) [H. Il5l|. 

The KSE is an indispensable tool in the modern de- 
scription of dynamical systems. Firstly, from it we learn 
about the degree of chaoticity because its inverse is the 
time scale of predictability. Secondly, this dynamical en- 
tropy is a well-defined quantity for both equilibrium and 
non-equilibrium systems. Thirdly, when tiny deviations 
of initial conditions that were not observable in the be- 
ginning are enlarged by the evolution, this can be inter- 
preted as the production of information about the ini- 
tial conditions. Finally, the KSE is known to be related 
to macroscopic properties such as transport coefficients 
[ll[ll[li,[li,[23,[2ll,[23. 

Our objective is to compute the KSE for the wet gran- 
ular gas. Pioneering work has been done by H. van 
Beijeren, J. R. Dorfman et al. [2^, [23| in the analytic 
treatment of sums of Lyapunov exponents for the gas of 
hard elastic spheres. We develop a generalization of the 
method suggested in [2^ . 

This article is organized as follows. In section |TT| we 
describe in detail the hysteretic interaction of wet gran- 
ulates. This Capillary Model allows the sticking of par- 
ticles by attractive forces in contrast to the "Standard 
Model" for dry granulates which assumes that a certain 
fraction of energy is lost instantaneously by inelastic col- 
lisions. In section IIIII we use the terminology developed 
in section|lT|to relate the behavior of the two-particle sys- 
tem to the full iV-particle system. Thereby we are lead to 
determine the probability distribution for colliding pairs 
of particles in section HVl In section |V| we derive the for- 
mula that expresses the expansion of velocity space as 
a function of the two-particle initial conditions for arbi- 
trary spatial dimension. In section IVII the results of the 
sections IIIUlV] are combined to accomplish the computa- 
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FIG. 1: Radial forces between a pair of wetted spheres. Solid 
line: The radial force of the Extended Capillary Model is 
plotted versus the center distance r. There is no interaction 
between the particles as they approach. After the collision ap- 



plies F(r) = -Fn, 



for r £ ((T, rcrit), otherwise there 



is no force. Dashed line: Experiments yield a decreasing force 
law [25I . [2^ with a discontinuity at the rupture. Therefore the 
even simpler Minimal Capillary Model which assumes a con- 
stant force that drops to zero at the critical separation is a 
good alternative approximation. The hysteretic interaction is 
the relevant property which is described by both the Minimal 
and the Extended Capillary Model. 



tion of the KSE. 



II. THE CAPILLARY MODEL 

There is an experimentally well confirmed Capillary 
Model for the dynamics of wet granulates, that will be 
applied here The system consists of hard spherical 
grains with equal diameter a and equal mass m. These 
are covered by a liquid film, so that every time two par- 
ticles touch, a liquid bridge is formed. The Capillary 
Model assumes that bridges are formed instantaneously. 
As we focus on the dilute gas, we may restrict our con- 
siderations to pair interactions. 

Experiments and computations [2^, yield a capil- 
lary force law that is excellently described by 



F = 



TTja cos ( 



1-t- 0.74 s -I- 1.25 s2 



(2) 



with the wetting angle ^w, the surface tension 7, and 
s = sy^a7T4)ridgc being the surface separation s expressed 
in the natural length unit ■\/Vbridgc/o' of the liquid bridge 
volume Vbridgo- The Capillary Model assumes that the 
bridge pinches off at a critical surface separation s = Scrit 
(i.e. at a distance rcrit = f+Scrit of the centers). To lead- 
ing order, the rupture distance Scrit equals the cubic root 
of the bridge volume Vbridgo- The energy that was stored 



in the stretched liquid bridge before the rupture is dis- 
sipated into the liquid and lost for the granular motion. 
We emphasize that this is the only dissipative mechanism 
in the Capillary Model (cf. the review article [s*], espe- 
cially Fig. 7 therein, for the capillary regime in which the 
Capillary Model applies.) In the moment of the rupture, 
the system is non-Hamiltonian because the atomic de- 
grees of freedom of the liquid to which energy flows are 
masked out in the description of the granular dynamics. 
Of course the forces acting on the grains are finite at the 
rupture, so that the trajectories (as functions of time) are 
continuous in the granular phase space and differentiable 
with respect to the initial state before the rupture. 

By a collision we denote the moment when two par- 
ticles in the entire TV-particle system touch each other. 
Since we are interested in statistical statements and a 
point in time is of measure zero, we can assume without 
loss of generality that there is a unique sequence of colli- 
sions. For a certain pair of colliding particles, we refer to 
the "collision cycle" as the time interval [ti,tf] that com- 
prises the collision of these two particles. The collision 
cycle starts at when the last particle of the two breaks 
free from its former collision partner and ends at t{ in the 
moment when the liquid bridge between them ruptures. 

During its collision cycle the radial motion of the two- 
particle system traverses a hysteresis loop. This is shown 
in Fig. [1] for the force ^ (dashed line) and for a simpler 
force law (solid line). The solid line in Fig. [1] falls off lin- 
early with the surface separation s. This is the Extended 
Capillary Model in contrast to the Minimal Capillary 
Model of Q which assumes a constant force. The cor- 
responding hysteretic "potential" of the Extended Cap- 
illary Model is 



-^loss 



-1, 

0, 

00, 



a < r before first collision, 
, a < r < Tcrit after collision, 

fcrit < r after collision, 
r < a. 



(3) 



In both, the Minimal and the Extended Capillary Model, 
the hysteretic loss of energy, i.e. the area -Eioss = 
-X, "Frdr in Fig.[Tl is a characteristic system prop- 
erty. When the energy in the center of mass system is 
below -Eiossi colliding particles will form a stable bound 
state with periodic collisions. With faster relative motion 
the liquid bridge exists for a finite time until the particles 
scatter off each other. We define a corresponding rela- 
tive velocity wioss by i?ioss = ™^^ioss/4 (with the additional 
factor 1/2 because rn/2 is the reduced mass). From this 
point on we distinguish between scattering events and 
collisions leading to bound states. For the scattering, the 
restitution coefficient e = Ef/Ei of the Capillary Model 
is an increasing function of the initial energy or velocity: 



1 



Ei 



or €{v[) 



loss 



(4) 



The binding threshold -Eioss of the Capillary Model con- 
trasts sharply with the widespread models for dry gran- 
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FIG. 2: The effective potential for Vi > vioaa and three dif- 
ferent impact parameters. For the solid line in the middle 
b and Vi are critical. For the higher b (fine dotted line) the 
particles are bound, for a lower b (roughly dotted line) they 
scatter. The inset shows the complete space of collision pa- 
rameters. The critical velocity v^it (plotted in units of vioss 
for rcrit = 2a) as a function of the scaled impact parameter 
6/(7 divides the plane in bound and scattering states. 



ules that assume either a constant or with increasing ve- 
locity decreasing coefficient of restitution for the collision 
of viscoelastic particles , [s^ . 

Let us denote by Ucrit the critical modulus of the rel- 
ative velocity Vi = vi — V2, that determines wether the 
incoming particles will form a bound state or scatter. For 
head-on collisions (impact parameter 6 = 0) Wcrit = floss, 
otherwise Wcrit > i^ioss since there is additional energy in 
the rotary motion. The next step is to determine Ucrit as 
a function of b. 



(of the liquid bridge potential given by ^) does not 
reach a maximum in r after the collision and leads to 
a bound state. For most Vi > v\oss the particles scatter, 
but there are some bound cases with high angular mo- 
menta, corresponding to high impact parameters. Fig- 
ure [2] shows three effective potentials for a given initial 
velocity Vi and different impact parameters b. In the 
case drawn with solid lines, b and Vi fulfill the critical 
relation Vi = Ucrit(6). For the higher b (fine dotted line 
in Fig. [2]) we have Vi < Wciit(6) so that a bound sys- 
tem is formed. Hence the criterion is that (jfcsir) touches 
the asymptotic energy -Eioss — mvf/4: in a single point. 
For the Extended Capillary Model it is possible to calcu- 
late these intersections explicitly. These are the roots of 
(i?ioss — mvf/4: + (j)cs) r^, which is a fourth order poly- 
nomial in r with one trivial root at r = and another 
unphysical root for r < a. So there are two real roots 
for the bound state which turn into a complex conjugated 
pair of roots for the scattering state. (Since the derivative 
of (pcS is continuous and negative at r = ^crit, the turn- 
ing point Tmax of a bound state follows correctly from 
this analytic consideration to be r„iax < ^crit without the 
need to take the non-analytic point r — r^it of 0off into 
account.) The easiest way is to compute the discriminant 
of the fourth order polynomial (E'ioss — mvf /4: + (pcs) r"^, 
which is equal to 



(8w^ - 4v'^ (57 + 9) + (27 + I87 - 7^)) ^ 

v'^ + 3v^-f + 3v^ (7 - 1) 7 + (7 - 3) 7^ - 7^ 



Determination of the Critical Velocity 

The bridge interaction is a central force problem. If Vi 
is lower than i;ioss, the effective potential 



4>cs{r) 



mb'^v? 
4^2 



(5) 



with 7 = fj 



The discriminant vanishes as the 



two physical roots coincide. Since the impact parameter 
b enters the problem only trough the angular momentum 
term in ([5]), the discriminant is a quadratic function of b-. 
Therefore it is elementary to give bcriti^i) as the inverse 
fimction of Vcnt(b) explicitly: 



b^nt{vi) \/-8 - 20(52 + (54 + 16w2 + 20d^w^ - 8u)4 - 5 (8 + ^2 _ gw^f^^ 



4:V2w{S - 1) 
I 



(6) 



with (5 = '''="* and w — . This function is plotted as 
inset in Fig. [2J Much more concise is the corresponding 



function for the Minimal Capillary Model: 

'^^loss 



Vcnt{b) 



1 



62 



(7) 
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FIG. 3: The collision sequence s{t) and the collision cycles: 
the step function s{t) is the total number of collisions in the 
entire A^-particle system until time t. The horizontal solid 
and dashed bars symbolize the collision cycles for scattering 
and bound pairs respectively. For the derivation is important 
that overlapping cycles affect different pairs of particles. The 
dashed arrow indicates a third particle that hits and breaks 
up a bound two-particle state. 



In the following sections including the main results ((4T|) - 
(H^ of this article, we shall be completely general without 
the need to specify for the Minimal or Extended Capillary 
Model. 



III. HOW TO RELATE THE TWO-PARTICLE 
SYSTEM TO THE iV-PARTICLE SYSTEM 

In the previous section we have shown how on the level 
of two-particle interactions the most important property 
of the real wet granular gas, namely the hysteretic bind- 
ing and breaking of liquid bridges, can be modelled. Fur- 
ther, we have seen that the bond energy of the liquid 
bridge gives rise to the sticking of particles. In this sec- 
tion we treat the many-particle system. 

Let V denote the mean collision frequency per particle. 
If the modulus of the initial relative velocity Wj is lower 
than Ucrit, so that particles stick together, the collision 
cycle is not terminated until a third particle bumps into 
the bound two-particle system. We assume that the out- 
state of such a three-particle event contains free particles, 
because the formation of higher mass clusters is rare in 
the gas-like state (cf. Fig. [TUl) . The pair interactions tak- 
ing place in the A'^-particle system may be envisaged as 
shown in Fig. [3l The number of collisions up to time t is 
denoted by s{t). Since s{t) is strictly monotonic its in- 
verse t{s) exists. The collision rate of the system, s/t{s), 
tends for s — > oo to Ni>/2 (each collision involves two par- 
ticles) . To have the steps visible Fig. [3] has been drawn 
for low N. The horizontal bars represent the concept of 
collision cycles introduced in the last section. There are 
two particles which are going to collide. As the begin- 
ning of the collision cycle we take the time when the last 
of these two particles has ruptured its liquid bridge con- 
nection to some previous collision partner. The collision 
cycle will end when these two particles rupture the liquid 



bridge between them. Thus a solid arrow in Fig. [3] shows 
that one of the particles which just finished its collision 
cycle immediately begins another one. The dashed arrow 
indicates that a third particle (that came out of another 
collision cycle) ends a bound two-particle state. 

With this picture in mind the computation of the KSE 
can be tackled. As stated by Pesin's theorem the KSE 
equals the sum of all positive Lyapunov exponents, be- 
cause the system is closed and sufficient chaotic • Lya- 
punov exponents describe the rate at which a certain di- 
rection in phase space grows or shrinks for large times. 
There is a orthogonal set of Lyapunov vectors describ- 
ing the direction while the associated Lyapunov exponent 
Xj describes the exponential rate 



O(t)^e,(0)e 



\i t 



(8) 



for long times t. According to the sign of Xj one speaks of 
stable or unstable directions. The deviations in the ini- 
tial conditions are infinitesimal small, i.e. the Lyapunov 
exponents characterize the tangent space map associated 
with a certain trajectory. In an ergodic system the Lya- 
punov spectrum {Xj} is independent of the trajectory 
according to Oseledec's theorem [23|. There is no doubt 
about the ergodicity of the gas of 3> 1 hard spheres 

Since in a dilute system the free flight time and the 
mean free path are large compared to the interaction 
time and the range rcrit of the interaction, perturbations 
of velocities are amplified as compared to spatial devia- 
tions [1^. This is not to be understood as a neglect of 
the spatial Lyapunov exponents. The Capillary Model 
is symplectic so that for each positive exponent Xj 
there is a negative exponent Xk = —Xj and the fact that 
the spatial deviations remain small means that the spa- 
tial directions mainly contain negative Lyapunov expo- 
nents, while the positive ones are assigned to velocities. 
So the conjecture is that the velocity space coincides (ap- 
proximately) with the unstable manifold of the system. 
Based on this conjecture the KSE, h^s, is given by the 
logarithmic volume growth rate in velocity space: 



^KS 



lim — — In 

s^oo t[S) 



det Y[ Mi 



(9) 



The deviation matrix Mi of the i's collision cycle is re- 
stricted to velocity space, so that it describes the evo- 
lution of velocity perturbations. There are three cru- 
cial points here: (i) This limit exists by virtue of Os- 
eledec's multiplicative ergodic theorem 27]. (ii) We have 
an unique collision sequence, (iii) Although there are 
pair interactions occurring with time overlaps, there is 
no ordering problem when writing down the total devi- 
ations as a product of collision cycles, because the co- 
existing liquid bridge interactions affect always disjoint 
pairs (by the assumption that there are two-particle clus- 
ters only) and deviation matrices of disjoint pairs com- 
mute. Therefore the matrices Mi can describe the full 
collision cycle of a single pair of particles, ignoring all 
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other interactions taking place simultaneously in the A'^- 
particle system. Our approach differs from ^ , because 
the Capillary Model has a hysteretic interaction with fi- 
nite interaction time. The dry limit follows by turning 
off the interaction, iSioss — > 0, as a special case. 
The expression (jH]) can be simplified dramatically: 



^KS 

N 



lim — !— In 



det Y[ Mi 



1 1 " 

- lim — - VlnldetMi 

N s^oo t{s) ^ ' 



N 
1 



1=1 

s E'-ilnldetM, 
— lim — — — 

TV s-^oo t[s) s 

^ (ln|detAf|) . 



(10) 



Herein the brackets < . . . > denote averaging over the 
two-particle phase space only. 

Since we expect the Lyapunov exponents to be of the 
order of the collision frequency i/, they are (according to 
the limit in ([9])) only well-defined if we let the system 
evolve for a time 



tL 



yapunov 



1 

^ — = icoU ■ 



In the subsequent discussion we will point out that this 
can be fulfilled even if there was no external driving 
mechanism to keep the dissipative system in a stationary 
state. Clearly, without a thermostat the system cools, 
T < 0, [iOjIiJ- The collision frequency v is of the order 
\T\/E\oss- On the other hand, cooling will be irrelevant 
on time scales below tcooi = T/\T\ . So the hierarchy 

^coll ^Lyapunov ^-cool 

of time scales can be fulfilled if 



(11) 



This implies that for weak liquid bridges as compared to 
the thermal energy we may speak of a Lyapunov spec- 
trum independently from the question of the thermostat. 
No additional limitation is set, since the condition pT|) is 
already required to be consistent with the gas state (dis- 
playing mainly single particles instead of clusters) which 
is studied in this work. 

Two tasks remain. The determination of the probabil- 
ity distribution for the formula (fTO|) is done in the next 
section. To make use of momentum conservation the sub- 
space is spanned by the center of mass position R = liila 

and velocity V = i^iiiS- of the two-particle system, as 
well as the distance r = ri — r2 between the centers of 
the spheres and their relative velocity v = vi — V2- The 
last step is to compute for any spatial dimension D the 
matrix M appearing in (jlOp . which maps for a specific 
point in the 4_D-dimensional phase space {R, f, V , v) the 



initial velocity deviations 



5Vi 

5v\ 



(12) 



from the beginning of the collision cycle to the final de- 
viations 



5Vf 
5v{ 



= M 



6Vi 

Svi 



(13) 



at the end of the collision cycle. This is done in section 

El 

Before we derive the joint probability density a com- 
ment on the velocity distribution itself is in order. It is 
well-known that for dissipative gases the velocity distri- 
bution can deviate from the Maxwell-Boltzmann veloc- 
ity distribution [29j depending on the state and driving 
mechanism. For explicit results we shall use the Maxwell- 
Boltzmann velocity distribution, 

(14) 



with 



2T 



. The result for the KSE will also be given 



in a form that is readily evaluated for any velocity distri- 
bution. For the distribution the modulus of the 
initial relative velocity is distributed according to 



r(f) ' 



dui 



(15) 



IV. THE ENSEMBLE AVERAGE 

We determine the probability distribution for two par- 
ticles under the condition that they will collide in the 
future. Therefore we depict the initial configuration 
of an arbitrary pair of particles in relative coordinates 
ri — ri— r2 as follows (Fig. [4|) : we rotate our coordinate 
frame such that the horizontal axis is per definition 



"Hi 

Vi 



(16) 



with the initial relative velocity Vi = vi —V2. This means 
that particle 2 rests in the origin while particle 1 moves 
horizontally to the right. Clearly, the particles will collide 
if and only if 



(i) the impact parameter is low enough, 
6= Wr,2 



n, 



< 



(a) and particle 1 is to the left of particle 2, 
in^Vi) < 0. 
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tide 2. 1 2 3 ^ 



For any pair of velocities vi,V2, there are initial relative 
spatial positions that lead to a collision. So we have to 
integrate over the entire velocity space x M^, 



FIG. 5: The distribution of Xi after averaging out the veloc- 
ities. The dashed curve is an exponential distribution with 
the same mean. Clearly P(xi) deviates from an exponential 
at distances Xi below the mean free path /. 



D 



d^'vi / d^w2 e 



(17) 



We take condition (i) into account by integrating the im- 
pact parameter over the interval [0, a]. From the conven- 
tional assumption of molecular chaos (i.e. the positions 
and velocities of two particles are uncorrelated) follows 
that the impact is uniformly distributed within the cross 
section, 



the probability density of the initial separation Xi is 

dxi f°° dx2 



P{Xi\vi,V2) = C 



I 



I 



-{xi+X2)/l 



X S { Xi — xi — ] S { — — — 



C e" 



"l+"2 



lD-2 



P{b) db^{D~l)- 



T-D-l 



, < 5 < cr 



(18) 



Further, we need to know the horizontal distance a;i > 
to the collision point. Together with the impact param- 
eter b this determines the relative spatial position com- 
pletely in the plane of incidence, since according to (m), 
r = b By — {x\-\- Vc^ — b"^) Cx always points to the left. 

The probability distribution of follows from the dis- 
tance covered by the particles in the laboratory frame. 
Denoting by xi and X2 the length that particle 1 and 2, 
respectively, have travelled in the laboratory frame since 
the beginning of the collision cycle, we have the equal 
time condition 



Xl_ _ _X2_ 

— f-frcc — 
Vi V2 



(19) 



where ifrce stands for the time of free flight that both par- 
ticles have in common. From this follows for the initial 
separation of particles 



Vl' 



-Xi . 



(20) 



The probability density of the travelled distances xi and 
X2 are known in a gas to be 



I 



J 



1,2. 



(21) 



The length scale I is the mean free path in the laboratory 
frame. Hence, under the assumption of molecular chaos 



up to a normalization factor. Obviously this yields the 
integration 



Vl + V2 



dxi 



(22) 



as part of the ensemble average. Putting pT)) . pB]) and 
together we can compute arbitrary expectation val- 
ues: 



(...) 



p-i)(-)7 d-.,/ d 



D V1+V2 
V2 



X e 



db dx^ 



(23) 



with Vi = 11^1 — W2||- In passing we take a look at the 
distribution of Xi in Fig. O The joint distribution ([^5]) 
implies that Xi is approximately distributed according 
to an exponential fall off, as one may expect, because 
the distances in the laboratory frame follow such a law. 
However there are differences: the mean is lower, e.g. 
< > « 0.71 I for D = 2, and the distribution falls off 
faster than exponentially for small Xi (cf. [23]). 



V. THE EXPANSION OF VELOCITY SPACE 

We aim to compute the determinant of the matrix AI 
as defined by Eq. There are always two distinct 

deviation matrices A/bound for vi < Wcrit and Mgcatt for 
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Vi > 'I'crit, SO that the phase space average naturally de- 
composes into 

(ln|detM|) = (ln|detMbound|)^,<„„,^ 
+ (ln|detM,eatt|)„,>,„,, . 

After determining these matrices, Eq. pO|) will enable us 
to compute 

^ = ^ [ (ln|detMbound|),,<„,,, 

+ (ln|detM,catt|).,>„,„J . (24) 

Because of momentum conservation, — Vf, the matrix 
M is of the blocked form 



M = 



Id Qd 



where Id and Qd are unity and zero matrices of dimen- 
sion D X D respectively. Therefore the only contribution 
to the growth in velocity space stems from the relative 
velocities, 



det M = det M' 
The final relative velocity [s^l is 



(25) 



Vf - 



- <ss (cosi9 + sinz? Cj,) . (26) 



As defined in (fT6|) ex points in the direction of the in- 
commg velocity and e„ = e^^ x ^'n = -rr^ — ^ i^^iu is 
the orthogonal vector spanning the plan of motion, such 
that 

r = -X-fix + bey 



with Xi = -)- Xcoi and ^coi = -(?\;oi, 4) = - 6^ is 
the x-distance of the particles in the moment of collision. 

When considering deviations of ((26|) one has to take 
into account contributions due to the change of the angle 
m i? = i9(6(r,7/i),^;), 



5b 



ddX- 

TIZ — """y 
db V, 



di9 



SVx 



(27) 



as well as contributions caused by rotations and inclina- 
tions of the orbital plane of motion: 



/ Sex \ 




6ey 




Sez 




\ ■■ ) 





I 



V 





SVz 





b 



b Vi 





Cj, 

ez 

.;[■■) 



(28) 

The Eqs. and ((281) hold for arbitrary spatial di- 

mension D. The resulting deviation matrix M' is rather 
complicated: 



M' 



( cop? _g^.^^gij^^ 

sif^ -hewii?„cosi? -f(l 




-Xii?6)esin7? 
Xii?fc)ecos7? 




e (cos I? -f sin i?) 



(cosi9 



X, 



11?) 



(29) 



with the restitution coefhcient ^ and the abbreviations 
i?fc = i?„ = The determinant of M (which equals 
M', cf. Eq. ([25]) ) is surprisingly simple: 



This reduces in the dry case, uioss = 0, to the expressions 
(18) (£"=2) and (19) (£'=3) in [H. The first factor in 



30|) is always non-zero since ||- < 



det M 



d_d_ 
db 



X ^1 + ^ sini9 
where we eliminated Xcoii ^ Xi using 

a;coii'?6 « -2 

a^coii . , „ 

— — sm « 2 — 2^ 

cos I? « 2^ - 1 . 



D-2 



(30) 



VI. RESULTS FOR THE KOLMOGOROV-SINAI 
ENTROPY 

In Fig. [6] the relative dynamic r(t) (which equals the 
motion of one of the two particles in the center of mass 
system up to a factor of 2) is sketched. In both cases, the 
determinant of M is of the form (jSO]) , but the meaning of 
the angle d{b, Vi) is quite different. For impact velocities 
above the critical value, i) is the scattering angle 



cib, Vi) 



arcsm ■ 



8 



(a) 



(b) 















/cnl 








/Vent 



















FIG. 6: The relative motion for (a) sticking and (b) scattering. 



:it 

whereas for Vi < ^critical the angle "i? is a function of time, 



TT . b 

arcsin — 

2 cr 



t3 



yarc(6, Vj) 
iarc(&, Vi) 



Here ^3 denotes the time during which the two-particle 
systems remains bound until it is freed by a third par- 
ticle. The angle between two contacts iparcib-,Vi) equals 
2 JJ"""''-'''^') dipif,{r) and there is a similar integral for the 
time tare it takes to run through one arc. The index 
(j) ought to remind us that the potential ([3]) enters only 
through these integral expressions. For ^3 ^ tare the 
angle inbound grows linearly with time, while the bound 
oscillations i^Sosc are negligible. 

Depending on the details of the interaction potential, 
iparc and tare Can grow beyond all bounds as the pair 
{b,Vi) approaches the critical line (&, ■ycrit(^)) (cf. Fig. (2) 
in the bound regime (from below). This singular behav- 
ior occurs in the Extended Capillary Model (linear force. 
Fig. [T]), whereas in the Minimal Capillary Model (con- 
stant force) both quantities remain finite. Close to the 
divergence the motion is an outward directed spiral, so 
that the turning point is never reached and the periodic 
collisions end. The interaction time can also diverge for 
scattering states (reaching the critical line in Fig. [2] from 
top) , but this singularity is integrable with respect to ve- 
locity. In the bound case the divergence is cut off by the 
third particle and because of angular momentum conser- 
vation we have the estimate 



■&honndit3,b, Vi) < const -|- 



hbvi 



(31) 



We will use the right-hand side as an approximation. 
The stopping time ts is a random variable itself and dis- 
tributed according to 



Vi 



dt3 



(32) 



for a given center of mass velocity Vi of the bound system. 
There is a smaller mean free path I' for the bound two- 
particle system: since its total cross section changes with 
time the effective diameter tJoff is |cr so that the mean 
center-center distance at contact is |ct. 

of 



Another factor 



I is caused by the mass ratio [31|, thus 



D-l 



(33) 



In the following, we shall evaluate averages that are linear 
in is, so that we can forthwith substitute the expectation 
value, — j^, of the distribution ([32]) . Then from ((3T|) 
follows 



db 



-ivi) 



Vi I' 
Vi (72 



(34) 



is at least of 



In both cases, binding and scattering 
the order of - , while Xi is of the order of the mean free 
path 



V2tt— 



(35) 



with n being the number density of grains. Formulas 
for the mean free path are well established [30] and other 
characteristic quantities for the motion of tracer particles 
are also available [Slj . We remark that investigating the 
trajectories of tracer particles is a promising technique 
for the experimental confirmation of results presented in 
this article. 

Our goal is to expand the KSE in the small dimen- 
sionlcss parameter na^ 1. So this is an expansion for 
the dilute wet granular system. The unity in the first 
and the last factor in pO)) contributes to the KSE only 
in linear and higher orders, while we are interested in the 
logarithmic and zeroth order terms: 



IdetAfl = 



db 



1 - 



■9{vi- Vcrit) 



/Xi_ 
V b 



sin'd 



D-2 



(36) 



With the step function 6, Eq. is valid for scattering 
and binding, because we assume that the collision with 
the third particle rethermalize the two-particle system, so 
that the next collision cycle starts with the same initial 
distribution. Since the 'third' particles have an energy of 
the order of the granular temperature T ^ i^ioss we can 
safely neglect the formation of bound states of three or 
more particles (cf. Fig. fTU)) . A cluster size expansion will 
be discussed at the end of this section. 

After introducing the appropriate length scales / and 
a we are lead to examine 



^KS 

TV 



{D - I) In 1 - (D - 2) /in - 
a \ 0" 
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(D-l) (In 



2^ 



In a 



db 
1 I Ine 



'Ui<l'crit 



D 

+ (L>-2)(ln|sini9|) 



In o- 



d-d. 



db 



fi>''crit 

(37) 



The first two terms in the square bracket yield 



Inna^ -Cd 



(38) 



InTT- 



with a numerical constant Cd 
^ - (L> - l)lnr (-^). This is independent of the 
ensemble average and the interaction potential. 

If Xi was distributed exponentially with mean Z, the 
third term in ([57)1 would give rise to the negative of Eu- 
ler's constant, — 7EuiGr ~ —0.5772, independent of the di- 
mensionality of the problem. As discussed before, lower 
values of xi are favored. That is why we find by numerical 
computation a lower expectation value, e.g. for D = 2: 



In 



Xi 



-1.01. 



(39) 



The fourth term in ^ is (cf. Eq. ([34| ) 

5l?bound 



In cr 



db 

Inna" + C 



1 



D 



l'i<1'crit 



l'i<'"crit 



/ l'i<Dcrit 



(40) 



with the numerical constant Cd = {D — l)lnj + ^ 



D-l 



2 ln^-lnr(i2±i: 



Together with (|38p the logarithm In na^ herein forms 
the leading term of the density expansion. Therefore the 
logarithm Inntr^ in ([40]) is a correction of the leading 
term as it is known for the dry case The KSE has 
the following density expansion: 



= -vADlnna" + uBD + 0{na°), (41) 



hKS 

N 

with the leading coefficient 

-^loss ^crit \ 



A, 



A, 



T 



<J J 



D-l D-l /m\f 



TO \ ■ 

r(f ) virJ 



dbb 



D-2 



rD-l 



t(b) 



D-l 



/Q O- JO 

and the density independent part 
Bd = i[-Cz3 + (i5-l)(lnf ) 



,(42) 



Cnd) 



+ (lnf 



(In (a 



I SlJs, 



db 

(In I sin I? I) 



) ) ui>l)crit 



(43) 
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FIG. 7: The increase AAd = ''g"'"' of the leading coefficient 
Ad ~ ^^-T^ + AAd- The solid line is for two, the dashed 



line for three dimensions D. Since A — 



in the absence 



of the liquid bridge interaction we recover the result for dry 
granulates as a special case. With the approximation for the 
wet granular gas used in the derivations one is restricted to 
temperatures above the bridge energy Eioss- Otherwise the 
method applied has to be extended to take clusters of more 
than two particles sticking together into account. The far 
extreme case, -Eioss ^ T, is known as the so-called sticky gas. 



The general form of the leading term, valid for any 
velocity distribution, is 



An = 



D- 1 



bound 



(44) 



We want to emphasize that so far all results of this 
section are general with respect to the spatial dimen- 
sionality of the problem and the details of the particle 
interaction. The probability Pbound = {^)vi<v„it ™ |44|) 
is given by integrating velocity and impact factor over the 
bound states in Fig. [21 Only here the detailed interaction 
models ([6]) and ([7]) enter the problem. 

Let us now turn to explicit results. For the Gaussian 
velocity distribution (jlSp and odd spatial dimensions the 
velocity integral of Pbound is an incomplete Gamma func- 
tion. In even dimensions the integral is elementary, yield- 
ing for D = 2 



^2(e,7) =1-5/0^^2; e- 



_ Bio 



as a function of the bridge energy over granular temper- 
ature, £, and the wetting content, 7 = rcrit/c > 1. The 
remaining integration variable is the impact parameter, 
X — h/ a. The excess of the critical energy over the bridge 
energy, f[x,^) = -Ecrit/^'ioss, depends on the model de- 
tails. In the Minimal Capillary Model from Eq. ([7|) fol- 
lows 



/(^'7)=(i-:^ 



The coefficient An oi the Minimal Capillary Model is 
plotted in Fig. [7] as a function of the liquid bridge energy 
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for two and three dimensions. Very similar curves follow 
from the Extended Capillary Model. For the plot the 
limit of short liquid bridges, Tcrit = f, was chosen. This 
corresponds to a small amount of liquid that is just suffi- 
cient to wet the surface roughness of realistic spheres. In- 
dependent of Tcrit/f > 1, in the dry limit (or equivalently 
the high temperature limit) Ajj approaches (D — l)/2, 
which is the known result for hard spheres [23l |. For a 
higher content of wetting liquid, rcrit/c > 1, the depen- 
dence of the leading term on the binding energy becomes 
flatter, but in an experimental situation there is a simul- 
taneous gain in -Eioss when liquid is added. Varying the 
surface tension of water by adding a salt to the wetting 
solution is an experimentally feasible way to measure this 
curve directly with a fixed amount of wetting liquid, such 
that rcrit/c can be kept constant. 

From this graph we see the sensitive dependence of 
the KSE on the cohesion force of the wetting liquid. To 
gain analytic insight we investigate exemplarily the two- 
dimensional case plotted. Substituting z = 1/(1 — x^) 
gives 



Bo 



A2{e,l) = 1 



dz 



(45) 



Splitting up the integration at z 
the non- analytic part. 



1 /e allows to separate 



^2(e,l) = l - 



•= dz 



(46) 



The first integral in (|46p can be expanded in powers of 
£ £ [0,1) since z > 1. The second integral equals 2 
for e ^ 0, while its first derivative has a logarithmic 
divergence: 



^2(£,1) 



1 



ln£ 



(47) 



The constant C is exp (-z)/4z-Mn2/2-}-(l-l/e)/4 « 
0.56. This shows that the slope of is vertical at i?ioss — 
0. 

Let us finally look at the next higher order term Bu 
of the density expansion. For simplicity we restrict our- 
selves to the case Z? = 2, so that 



In^ 
Vi 



In a 



OA. 



db 



The last term in 
of dry granulates. 



"crit 

(48) 

■Ui>1'crit 

is exactly equal to unity in the limit 



lim ( In I (7 



db 




^loss 



FIG. 8: The coefficient B2 of the density expansion (|4ip . 
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FIG. 9: The two-dimensional KSE as a function of the density 
for three different bridge energies -Eioss. This energy depends 
on the amount of wetting hquid added to the granular gas as 
is indicated in the plot. Another way to change Bioaa is to 
add a salt or a surfactant. 



(T 



1-^ 



but decreases as the critical velocity increases when we 
turn on the liquid bridge interaction. The coefficient B2 
for the zeroth order in the expansion (14111 is plotted in 
Fig. [H It is known for the dry limit [23| , that the ac- 
cordance oi B]j with numerical simulation cannot keep 
up with the successful confirmation oi A^- The origin 
of this discrepancy is the assumption that the unstable 
manifold coincides with velocity space and it is quite in- 
volved to improve on that [32l |. In the dry limit our 
method yields B2 = —0.52(8), which is lower than the 
analytical estimate (S2 = 0.1045) and the simulated re- 
sult {B2 = 0.679) of m. 

From the knowledge of the coefficients Ad and Bjj 
follows the KSE in the dilute system for various wetting 
contents as shown in Fig. [Slfor D = 2. 
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The summation can be written as a systematic expansion 
in the cluster size: 



12 3 

Liquid Bridges per Particle 

FIG. 10: The probabihty for a sphere to have a certain num- 
ber of liquid bonds ending on its surface. This distribution is 
measured in a three dimensional molecular dynamics simula- 
tion of a wet granular gas with an occupied volume fraction 
of 3.9%, which corresponds to ncr^ = 0.074. The granular 
temperature T has been varied as indicated. The probabil- 
ity for two liquid bridges ending on one particle, as necessary 
for a three-particle-cluster, is suppressed by more than three 
orders of magnitude. An analytic approach to the KSE is fa- 
vorable because the direct numerical integration suffers from 
high computing times for the full tangent space dynamics and 
yields noisy results [33l . Issj . The liquid bond distribution 
shown is a robust and reliable single-particle quantity. 
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with the events Ti and T2 considered before in (|49|1 . 
The events Tj with j > 2 result in new many-particle- 
clusters which are exponentially rare components of the 
wet granular gas as is evident from Fig. [TOl We re- 
mark that the scattering of a bound state (T5) pro- 
longs the mean bond time t^, to become t'^ — at^, with 
a = 1 + 2Pt, + 3P|^ + ... = 1/(1 - Pts)^. The 
unity in front of this series corresponds to breaking the 
bound state in its first collision (T2), the second term 
corresponds to one scattering event of the bound pair 
and the following terms to multiscattering. The contri- 
bution to the KSE is proportional to the logarithm of 
this time, Intg = ln<3 — 21n(l — Pts). The first term 
Inia cx — ln{na^) is the wet granular contribution to the 
leading coefficient A as identified in Eq. (44). The sec- 
ond term gives a correction to the _B-coefficient which is 

of the order Pts = O ( \/E\oss/T 1 for three dimensions. 



The Cluster Expansion 

In Eq. (|24|) we considered events including bound 
states of two particles (a-|-b-|-c ab + c — > 

a + b + c) and scattering events (a + b a + b) by 
writing 

(ln|detM|) = (ln|detMbound|>„,<„„,^ 

-f(ln|detM,eatt|)l,>ri, ■ (49) 

The first term is proportional to Pbound which led to 
Eq. Here we wish to point out how to generalize 

the computation of the KSE to include clusters of higher 
particle number. All equalities in (10) hold for arbitrary 
types of events, when Mi denotes the deviation matrix 
associated with the ith event and v is the generalized 
event frequency. Referring to the event type by T we 
reorder the averaging. Collecting the events of type T by 
introducing (5typc(j),T (which is unity for an event T and 
otherwise zero) we write (. . .)rp for ^. . . \ype{j),T)- 




= (ln|dctM|) = ^ (ln|detMT|)T • (50) 



VII. CONCLUSIONS 
Summary 

We worked out the crucial difference in the interac- 
tion of wet granulates compared to the dry case. There 
is a liquid bridge causing a radial hysteretic force over 
finite distance. The detailed distance dependence is of 
minor importance. The decisive ingredient in the Cap- 
illary Model is the extraction of a bridge energy that 
is independent of the initial velocity in contrast to the 
"Standard Model" using a restitution to extract a cer- 
tain fraction of energy. 

We found an enhanced chaotic behavior of the wet 
granular system. The leading term in the expansion of 
the KSE with respect to the small density (ncr^ ^1) 
changed due to the possible sticking of particles. One 
can think of the prolonged interaction time enforcing the 
exponential separation in velocity space. The continuous 
but in general not differentiable transition to the limiting 
dry case has been established. 

This dynamical property recommends the wet granular 
system as a suitable candidate for experimental, numeri- 
cal and analytic tests of the Gallavotti-Cohen fluctuation 
theorem j34| which requires hard chaos. 
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Outlook 

In this analytic work we used an assumption on the 
unstable manifold and we neglected correlation effects 
in consecutive collisions. Although physically motivated, 
the next challenge will be to verify these assumptions by 
direct numerical simulations. 

The rigorous derivation of phenomenological laws such 
as the Navier-Stokes equation for viscous flow and the 
Fourier law for heat transport is a fundamental prob- 
lem under intense discussion. Relations between the Lya- 
punov spectrum of the microscopic dynamics and macro- 
scopic properties such as viscosity and heat conductivity 
have been established within the last years, most detailed 
for the Lorentz gas [H, [13, [M llS!, (20i l2ll (2^ . The sever- 
ity and importance of these relations become apparent 
from the fact they have to bridge the gap between micro- 
scopic reversibility and macroscopic irreversibility chal- 
lenging physicists since Ludwig Boltzmann. 

The dynamics of the wet granular system studied in 
this work follows a mesoscopic law including dissipation, 



and kinetic theory has already been extended to dry gran- 
ular matter The next step is to extend also these 
transport relations. We hope that our results on the Lya- 
punov exponents might stimulate this development. On 
the experimental side mechanical properties of wet gran- 
ulates are presently under investigation 

A further interesting problem is the computation of 
the KSE for dense wet granulates. This might lead to 
a novel description of clustering - as a non-equilibrium 
phase transition - in terms of the Lyapunov spectrum. 
Yet this problem is challenging as it needs new concepts, 
because the identification of the velocity space with the 
instable manifold is limited to the dilute gas. 
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